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S Abstract 

A continuum model for self-organized dynamics is numerically investigated. The 
j>. model describes systems of particles subject to alignment interaction and short- 

range repulsion. It consists of a non-conservative hyperbolic system for the density 
and velocity orientation. Short-range repulsion is included through a singular pres- 
(f-) sure which becomes infinite at the jamming density. The singular limit of infinite 



pressure stiffness leads to phase transitions from compressible to incompressible dy- 



namics. The paper proposes an Asymptotic-Preserving scheme which takes care of 
the singular pressure while preventing the breakdown of the CFL stability condition 
near congestion. It relies on a relaxation approximation of the system and an elliptic 
formulation of the pressure equation. Numerical simulations of impinging clusters 
^ show the efficiency of the scheme to treat congestions. A two-fluid variant of the 

model provides a model of path formation in crowds. 
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1 Introduction 



The modeling of self-organized dynamics is at the core of an intense scientific activity. 
Self-organized dynamics occurs in systems of active agents such as bacterial suspensions 
and sperm [SB] , insect swarms [12] , fish schools [U, [Tj5] , bird flocks [U HI] , mammal herds 
or pedestrian crowds J2U E7J HI] . These systems exhibit large-scale coherent structures 
that spontaneously emerge from the interactions between the agents but are not directly 
encoded in the interaction rules. Many mathematical models have been proposed to 
account for the emergence of large-scale order [H [9l [13l [181 UHl [301 ESI S21 HH] - One 
of the most paradigmatic models is the Vicsek particle system [51] which consists of self- 
propelled particles interacting through alignment interaction. Each particle moves with a 
constant speed and updates its direction so as to align with its neighbors up to a certain 
noise level. The Vicsek model exhibits phase transitions from disordered to fully aligned 
states when the noise is decreased or when the density is increased. 

In practice, more refined flocking models must be used to reproduce flocking behavior 
observed in nature. Among these refined models, the three zone model [TJ [Tj5] offers a 
good compromise between physical accuracy and simplicity. It considers three kinds of 
interactions: long-range attraction, medium-range alignment 'a la Vicsek' and short-range 
repulsion. In particular, short-range repulsion is believed to play a key role in the observed 
morphogenetic features of self-organized dynamics. The flock can indeed be seen as the 
region where all the particles are as close as they can be, given the short range repulsion, 
i.e. as the region where the particle density reaches a maximum value. The existence of 
a maximal allowed value of the density resulting from volume exclusion is referred to as 
the congestion constraint. 

Macroscopic models of fluid type (or hydrodynamic models) are particularly well suited 
to the modeling of large-scale self-organized structures by contrast to particle models (also 
known as 'Individual-Based Models' or IBM) which focus on the inter-particle interaction 
scale. Indeed, the numerical complexity of IBM's increases with the number of individ- 
uals faster than linearly. They become computationally very intensive for large particle 
systems which are the most relevant ones for the observation of self-organization. Addi- 
tionally, in some cases, self-organization can be directly encoded into the fluid model and 
facilitates the observation of the resulting morphogenetic features. In particular, this can 
be realized with the congestion constraint which results in a clear-cut phase transition 
between compressible and incompressible fluid regions, the latter being those where the 
congestion constraint is reached. 

The goal of the present work is to study a hydrodynamic model of self-organized 
dynamics with a congestion constraint. The starting point is the hydrodynamic model of 
[25] . which will be referred to as 'Self-Organized Hydrodynamics' or in short, SOH. This 
model is derived through a macroscopic limit of the Vicsek particle system [51]. The SOH 
model resembles the isentropic compressible Euler equations with two major differences. 
First, the local mean velocity f2 is constrained to have unit norm |f2| = 1 (specifically, f2 
is the mean velocity direction rather the mean velocity itself), and second, the pressure 
gradient term in the momentum balance equation is projected onto the plane normal to 
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f2, in order to preserve the norm constraint. This makes the system non- conservative and 
non-classical since the state variables belong to a curved manifold. 

The SOH model only takes into account the alignment interaction. In the present work, 
we take into account the short-range repulsion by considering a singular pressure when the 
density p approaches the congestion density p* . The SOH model with volume exclusion 
constraint has been theoretically investigated in [27J where the limit of an increasingly 
steeper singular pressure has been studied. This limit process results in a model where 
clear-cut phase transitions between uncongested and congested areas occur, the former 
corresponding to a compressible model and the latter, to an incompressible one. In [27] . 
the dynamics of the boundary between these areas (or boundary of the flock) has been 
studied. 

The present work is devoted to the derivation and study of adequate numerical meth- 
ods for this model. In [22], a numerical method for the classical isentropic compressible 
Euler equations with congestion constraint has been derived. The key point is the use of 
an 'Asymptotic-Preserving' (AP) scheme previously developed for the small Mach-number 
limit of compressible flows [T71 12"5| W7\ . Here, we adapt this scheme to the constrained 
SOH model. This adaptation is not straightforward, because of the non-conservative 
character of the SOH model and of the unit norm constraint on the velocity. In [33], it 
has been shown that a relaxation method is best suited to the numerical resolution of the 
unconstrained SOH model. The present paper shows that the relaxation method of [33] 
and the AP methodology of [22] can be combined in order to resolve the constrained SOH 
model. 

As an application of the presented methodology, a model of path formation in crowds 
is presented. The spontaneous emergence of trails in pedestrian counterflows in corridors 
is a well documented phenomenon [3JJ 133]- It has inspired several kinds of modeling 
approaches, both at the discrete [331 EH] and continuous [21 H2] levels. Our model intends 
to describe path formation in a crowd of people heading towards opposite directions. 
It requires a two-fluid extension of the model to account for the existence of two kinds 
of pedestrians having opposite target directions. Our model shows that the congestion 
constraint prevents the flow to cross through high density regions, forcing the emergence 
of paths of oppositely moving fluids. 

Throughout the paper, the key property of the developed scheme is the Asymptotic- 
Preserving (AP) property. It is defined as follows. Consider a continuous physical model 
M. £ which involves a perturbation parameter e (here, e is the parameter describing the 
stiffness of the pressure and Ai £ represents the SOH with stiffened pressure) which can 
range from e = 0(1) to e < 1 values. Let Ai° be the limit of Ai e when e — > (here 
M° is the constrained SOH model). Let now Ai% be a numerical scheme which provides 
a consistent discretization of Ai e with discrete time and space steps (At, Ax) = A. 
The scheme M.% is said to be Asymptotic-Preserving (AP) if its stability condition is 
independent of e and if its limit Ai% as e — > provides a consistent discretization of the 
continuous limit model The AP property is illustrated by the commutative diagram 
of fig. [T| The literature about AP schemes is recent, yet increasingly abundant and 
applied to various contexts (see e.g. [131 1321 [36J ) . 
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e -> 



A4 e 



A4° 



A -> 



A -> 



-Ma 



E -> 



-Ma 



Figure 1: Asymptotic-Preserving (AP) property: the upper horizontal arrow translates 
the assumption that the continuous model Ai £ tends to the limit model Ai° when e — > 0. 
The left vertical arrow expresses that Ai% is a consistent discretization of Ai £ when the 
discretization parameter A — > 0. The lower horizontal arrow indicates that the scheme 
Ai% has a limit M.\ when e — > for fixed A. Finally, the right vertical arrow expresses 
the AP-property: it says that the limit scheme A4° A is a consistent discretization of the 
limit model Ai° when A — > 0. 



Constrained hydrodynamic models obtained as limits of hydrodynamic models involv- 
ing a singular non-linear pressure have first been proposed in [TT] and studied in [7] for the 
compressible Euler equations. They have then been extended to traffic flow to describe 
the dynamics of traffic jams [8j. The SOH model is derived from a kinetic formulation of 
the Vicsek model (25]. This kinetic description has been shown to be a valid description 
of the particle system in [6]. Variants of the SOH model taking into account non-isotropic 
observation [31] . phase transitions [211 E2], viscosity [211122], trajectory control by curva- 
ture [2S] and precession [23] can be found. Alternate hydrodynamic models derived from 
the Vicsek system can be found in [20l 46J as well as related models in [151 IT6] . 

The paper is organized as follows. The model framework is given in section [2j Then, 
the numerical method is introduced in section [3j Numerical results are given in section 
|4j The application to path formation in crowds is reported in section [5j A conclusion is 
drawn in section [6} Finally, for the reader's convenience, the expression of the scheme in 
the two-dimensional framework is exposed in Appendix 1. 

2 Model framework 

The Self-Organized Hydrodynamics (SOH) model is written as follows: 



where p(x,t) e R + denotes the mass density, Q(x,t) 6 IR 2 , the mean velocity and 
q(x,t) G M. 2 the momentum density, depending on the spatial position a; £ I 2 and the 
time t > 0. The constants are such that A > and cel. The term p £ (p) = sp{p) is the 



Pt + Va; • q = 0, q = pf2, 

P (n t + c(n ■ v x )ci) + x (id - n g) n) v* {p £ (p)) = o, 
|n| = i, 



(2.1) 
(2.2) 
(2.3) 
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pressure, with a function p(p) specified below. We denote by ® the tensor product of two 
vectors. The matrix (Id — f2 <g> ft) is the orthogonal projection onto the plane orthogonal 
to fl. We do not specify any boundary conditions for the time being and suppose that 
initial conditions po and f2o are provided such that |O | = 1 and po > 0. 

Eq. (2.1) is the mass conservation (or continuity) equation while (2.2) is the momen- 
tum balance equation. We note that the model is similar to the isentropic compressible 
Euler model of gas dynamics but for two important differences. First, the velocity f2 is 
constrained to have unit norm, i.e. f2 e S 1 , where S 1 is the one- dimensional sphere. As 
a consequence, fit + c(f2 ■ V^fi), which is a derivative of f2, is orthogonal to f2. This 
is why the pressure gradient term AV x (p £ (p)) is multiplied by the projection operator 
(Id — ft Cg) ft) . In this way, the constraint (2.3) is satisfied at all times provided it is 
satisfied at initial time. This constraint expresses that f2 should be seen as the veloc- 
ity direction, or polarization vector (see [3J for another example of an active fluid model 
written in terms of its polarization vector). The actual velocity is a constant times the 
velocity direction. After some time rescaling, this constant can be chosen equal to unity 
(see |25j for details). 

The second difference is the presence of the constant c in front of the velocity transport 
term (Q-'V x )fi. In the isentropic Euler equations, c = 1, i.e. the coefficient of the velocity 
transport term is the same as that of the mass transport term V x ■ (pfi). In the SOH 
model, these coefficients are different, and may even be of different sign [31]. The difference 
reflects that, in swarming systems, velocity is not simply transported along with the fluid 
velocity itself. If c < 1, velocity propagates upstream the fluid. This situation occurs in 
car traffic, where drivers update their velocity according to the state of the flow in front 
of them. One may also find c > 1, which occurs when particles update their velocity 
according to the flow behind them [51] . This occurs e.g. in locust swarms because locusts 
want to prevent themselves from being bitten by their congeners behind them. These 
considerations illustrate the loss of Galilean invariance of the model when c ^ 1 and the 
consecutive anisotropy of wave propagation described by these models. A discussion of 
the lack of Galilean invariance in flocking models and its consequences on the propagation 
of sound waves can be found in [3D] . 

In the SOH model derived in [25], the pressure is just a linear function of the density. 
However, short-range repulsion can be included by introducing a non-linear pressure which 
blows up at the congestion density p*. This non-linear pressure relation has been derived 
from a microscopic dynamics involving finite-sized particles undergoing short-range repul- 
sion in [27]. The congestion density p* corresponds to the jammed state where particles 
are in contact to each other. Additionally, we assume that the intensity of the repulsion 
force is very small unless the density p is very close to p*. The parameter e characterizes 
the width of the density range where the pressure is finite: p £ (p) is supposed almost zero 
except in small neighborhood of p* and becomes infinite at exactly p = p*. This sudden 
'turn-on' of the repulsion force results in a sharp phase transition between uncongested 
and congested regions as shown below. 
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These modeling hypotheses are translated in the following way: p 6 (p) is defined by 

p £ (p)=ep(p), (2.4) 

where p(p) is such that p(p) — > oo, when p — > p*. To be more specific, we consider: 

1 



P(p) 



1 i_ 

p p* 



7 >0, 



(2.5) 



i.e. p behaves like usual isentropic fluid pressure at low density: p(p) ~ p 7 when p — > 0. 

<C p* while p £ (p) = 0(1) for p clo; 
displays typical functions p and p £ 



From (2.4), (2.5), it follows that p £ (p) = 0(e) f or p <C p* while p e (p) = 0(1) for p close 
to p* (more precisely for p = p* — (^(e: 1 / 7 )). Fig. 
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(a) p(p) and p e (p) (b) p £ (p) and pf (p) 

Figure 2: Left: the pressure p and the rescaled pressure p e as functions of the density p. 



Right: the two components p%(p) and p\(p) as functions of the density p (see section 3.2). 



Solutions of (2.1)-(2.3) depend on e and are denoted by (p £ ,Q £ , q e ). Their formal 



limits as e — >■ are supposed to exist and are denoted by (p, f2, q): 

(p £ ,n £ ,q £ )^(p,n,q) as £^0. 

When e — > 0, the singular pressure term p £ (p £ ) tends to different values according to 
whether the limit satisfies p < p* (uncongested region) or p = p* (congested region). We 
denote by C(i) and U(t) the congested and uncongested regions respectively and by I(t) 
the interface between these two regions: 

U(t) = {xeR 2 , p(x,t) < p*}, 
C(t) = {xeR 2 , p(x,t) = p*}, 
l(t) = dU(t) = dC(t). 

Uncongested region: We consider a point x G U(t), i.e. such that p(x,t) < p*. 
In this case, p(p £ ) —> p(p) < oo and consequently p £ (p £ ) = sp(p £ ) — > 0. Therefore, the 
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formal limit is the presureless gas dynamics with norm-one constraint on the velocity: in 
U(t), (p,Q,q) satisfies: 

Pt + ■ q = 0, q = pfl, (2.6) 

p (n t + c (n ■ v^n) = o, (2.7) 
|n| = 1. (2.8) 

Congested region: We now consider a point x G C(t), i.e. such that p = p*. In this 
case, p(p £ ) — > 00 and lim e _ s . (^p(p £ )) is undetermined. We assume that 

lim ep(p £ ) = p with < p < oo almost everywhere. (2-9) 

Now, p becomes an additional unknown of the limit system: in C(t), (p, Q, q,p) satisfies: 

p = p\ q = p*n, (2.10) 
V x -fl = 0, (2.11) 

p* (n t + c (n ■ Va.)n) + a (w - n ® n) v x p = o, (2.12) 

= 1. (2.13) 

The vector field f2 is divergence-free and of unit-norm. Smooth vector fields like f2 have 
a remarkable geometry (see Fig. |3(a)| and [27]): curves normal to fl are straight lines 
and f2 is constant along these lines. We denote by M the family of straight lines normal 
to Q. The quantity p is now an unknown like the pressure in standard incompressible 



models. The governing equation for p is obtained by taking the divergence of (2.12) and 



using (2.11) to eliminate the time derivative. It is written for x G C(t) : 

- AV* ■ ((Id -ft® ft) V x p) = cp*V x ((ft ■ V x )fl) . (2.14) 

This leads to a family of one- dimensional elliptic equations for p posed along the lines of 
the family M. Their inversion requires the knowledge of the boundary values of p on Z(t). 
These boundary values, as well as the velocity of the interface X{t) have been determined 
in [27]. They are linked to the jumps of p and f2 across X{t) and are found by solving a 
Riemann problem in the direction normal to I(t), supposing Z(i) is smooth. However, if 
X(t) is not smooth, which happens for instance at times when its topology changes, these 
quantities cannot be analytically known from [27]. Since we target a general situation 
where topology changes can occur, we will not use the conditions of [27] and do not need 
to recall their expression. 



The SOH model (2.1)-(2.3) is hyperbolic, with characteristic speeds given by: 



f ± = cos 9 + i ^(c - 1) cos 9 ± yj(c - l) 2 cos 2 9 + 4Xep'(p £ ) sin 2 , (2.15) 



where 9 is the angle between f2 and the propagation direction. If p e — > p* such that (2.9) 



holds, then, using (2.5), we get 

ep'{p £ ) = £>0r 1/7 ) -> oo. (2.16) 
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Lines of the family A/*(t) 




z(t) 

(a) Congestion region 




(b) Singularity of Q 



Figure 3: The congestion region C(t) at a given time t is enclosed by the curve X(t). Left: 
the vector field Q(x, t) is constant along lines normal to itself. These lines form the family 
M(t). The quantity p(x,t) satisfies one-dimensional elliptic equations posed along these 
lines. It requires boundary conditions along the interface X(i). Right: a situation where 
the vector field Q(x,t) has singularities 



Therefore, the characteristic speeds (2.15) tend to ±00. This explains why the limit 



problem (2.11 )-( 2.12 ) has elliptic characteristics, reflected in the elliptic equation (2.14) 



for p. The second term of (2.15) is the sound speed 

c s = i f(c - 1) cos# ± yj(c - l) 2 cos 2 9 + AXep' (p £ ) sin 2 9 

while the first term is the fluid velocity projected on the propagation direction u = cos 9. 
We can define the Mach number as the ratio of these two velocities 



M 



u 



cos 9 



°s | ( (c - 1) cos 9 ± ^(c - l) 2 cos 2 + 4A£p'(p E ) sin 2 9 



Due to (2.16), we generically have M. — > when e — )■ 0. Thus, in the congested region, 



the £ — > limit is a small Mach-number limit. We note that in the SOH model, the sound 
speed is anisotropic and depends on the angle between Q and the propagation direction. 
This has already been noticed in [50] for other hydrodynamic models of flocks. 

The uncongested and congested models can be unified by writing the limit system in 
the following way: for i6l 2 , 



Pt + ■ q = 0, q = pfi, 

p (sit + c(n • v x )n) + a (w - n 

\Q\ = 1, 
(p-p*)p = 0. 



n) = o, 



(2.17) 
(2.18) 
(2.19) 
(2.20) 



However, this expression does not specify what happens when the constraint (2.20) shifts 
from p = (in U(i)) to p = p* (in C(t)) at the interface X(£). Therefore, it must be 
complemented by interface conditions like those of 
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Remark 2.1 The pressureless gas dynamics model is not strictly hyperbolic. This results 
in some unpleasant features such as the formation of singularities or the generation of 
weak instabilities [TUf . According to the physical context, it is possible that a non-zero 
pressure prevails in the uncongested region U. This is easily modeled by considering a 
pressure p e (p) of the form 

p £ (p)=ep(p)+p B (p), (2.21) 



where p is given by (2.5) as previously and the background pressure Pb(p) is independent 
of e and has finite limit Pb{p*) when p — )■ p* . For instance, we can take 



Pb{p) = ^P 1 



(2.22) 



with the parameter k controlling the strength of the background pressure. In this case, the 



limit problem (2.6)-(2.8) in the uncongested regionU is changed to the following form : 



Pt + ^ x - q = 0, q = pft, 

p(n t + c(n-v x )n) + \{id- 
ini = i, 



n®n)v x p B {p) = o, 



which is nothing but the SOH model with pressure p B - In the congested region C, there is 



Pb(p*) < p < oo instead of (2.9). 



no change to system (2.10)-(2.13) except that now the pressure p satisfies the inequality 



The numerical resolution of problem (2.17)-(2.20) is quite challenging since the inter- 



face T(t) is a moving free boundary. It requires two different fluid solvers: a compressible 
one in the congested region C(t) and an incompressible one in the uncongested region U(t). 
The numerical discretization of the interface T(t) must be matched with the mesh used 
inside each of the domains C(t) and U(t). Since T(t) is a moving interface, some remeshing 
must be performed at each time step. Alternate methodologies such as level-sets may be 
envisionned. In all cases, these methods are computationally intensive and necessitate 
fine parameter tuning. Additionally, interface motion is induced by fluid motion, which 
itself is determined by the boundary conditions for p along the interface. However, in 
the case where T{t) changes topology, the formal asymptotic procedure of [27] does not 
provide any information about these boundary conditions. In this circumstance, it is not 
known what interface motion takes place. Then, there is no other possibility than solving 



the original perturbation problem (2.1)-(2.3), with a very small value of e. This is the 



strategy developed in the present work. 



To solve the perturbation problem (2.1)-(2.3) with an arbitrarily small value of e, an 



Asymptotic-Preserving (AP) scheme is necessary. Indeed, the limit e — >• corresponds to 
a small Mach-number limit. With a classical scheme (such as a shock-capturing scheme), 
the CFL stability condition leads to a time-step constraint At < C £ Ax with C £ — > 
as e — > 0. Therefore, a classical scheme is impracticable for small values of e. The AP 
property recalled in Figure [T] guarantees that, for small values of e, the scheme provides 
a consistent solution to the limit system (2.17)-(2.20). 
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In this work, we adopt the AP-method of [TTJ [2H1 HZ] which has already been used for 
the Euler system with congestion in [22]. However, the SOH model is a non-conservative 
model. In general, non-conservative models may have multiple weak solutions of the shock- 
wave type. Additional information is necessary to single out the correct shocks: a suitable 
path in state space connecting the two states on each side of the shock must be defined 
[3"§] . The choice of the correct path depends on some microscopic information which is not 
included in the macroscopic model and which, in general, is not available. Here, we look 
for solutions of the SOH model which provide good approximations of the underlying 
microscopic Vicsek model [25J. Such solutions are considered as being the physically 
consistant solutions. In jl3], it is shown that these solutions are well-captured by a 
relaxation model. It consists of a gas-dynamics type conservative system for the density 
and momentum supplemented by a relaxation right-hand side which drives the velocity 
towards a vector of unit norm. In the present work, we rely on this relaxation model as 
the one which produces the physically consistant solution and propose a modification of 
the method of jl3] which makes it AP in the small Mach-number limit. This modification 
is based on [2"2l |2"%] . This numerical method is presented in the next section. 

3 The numerical method 

3.1 Relaxation method and splitting 

We first present the relaxation method of jl3]. We consider the following relaxation 
system: 

d t p^ + V x -q £ >P = 0, (3.1) 
+ cV x ■ ( g£ y ) + AV. (p% P ^)) = i(l - 1^1 V'' 3 , (3-2) 



where the geometric constraint |f2| = 1 is removed. System (3.1), (3.2) is a conservative 
model with a relaxation term. We have: 



Proposition 3.1 j^j The relaxation model (3.1) -(3. 2) formally converges to the SOH 



model (2.1)-(2.3) as (3 -> 0. 



The left-hand side of (3.1), (3.2) coincides with the isentropic Euler system if and only if 



c = 1. If c 7^ 1, the characteristic speeds are given by 




c*-c) q - 2 +\E P '{p). (3.3) 
P 2 



where q is the component of the momentum along the propagation direction. These 
characteristic speeds are real and distinct if c > 1 (provided that p > 0) and may be 
complex if c < 1 for small values of sp'(p). Therefore, although the limit system (the 
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SOH model) is always hyperbolic, the relaxation system (3.1), (3.2) is strictly hyperbolic 



only for c > 1. We will see that, in spite of being non-hyperbolic, the relaxation system 
leads to a valid numerical solution of the SOH model in the relaxation limit when c < 1. 

Remark 3.1 In fj^j , o, linear pressure was used and the relaxation system could stay 
hyperbolic even with c < 1 for suitable choices of the parameter A. Here, the non-linear 
pressure and the smallness of the parameter e require c > 1 to guarantee hyperbolicity of 
the relaxation system for all values of p > 0. 



The plan is therefore to discretize the relaxation system (3.1), (3.2) with arbitrarily 



small values of (3. To this aim, we use the splitting method of [43]. Let At be the time 
step and (p n , fl n , q n = p n fl n ) be an approximation of (p, fl, q) at time t n = nAt. To pass 
from (p n , fl n , q n ) to (p n+1 , fl n+1 , q n+1 ), we first solve the following conservative equation 
with singular pressure p e : 



d t q + cV 3 



dtp + V x ■ q = 
q <E> q 



P 



+ AV a ( P e (p)) = 0, 



(3.4) 
(3.5) 



over the time-step At with the AP-scheme developed in [22] and initial conditions (p n , q r 
This leads to intermediate values (p n+1//2 , q n+1 ^ 2 ). Then, the relaxation system: 



Pt = 0, 



q 



(3.6) 
(3.7) 



is solved over the time-step At with initial condition (p n+1 / 2 , q n+l / 2 ) and yields (p n+1 , q n+1 ) 
Finally, we set f2 



n+l 



-,n+l I pn+1 



3.2 Time semi-discrete scheme 



We first focus on the time semi-discretization of the conservative step (3.4), (3.5). Fol- 
lowing [22], it is written as follows: 



n+l/2 _ n 

9 Af 9 + V, • q^l* = 0, 

gn+l/2 _ qn 



At 



+ cV a 



q n ®q n 



AV,(pg(p n )) + AV,(pKp n+1/2 )) = 0. 



(3.8) 
(3.9) 



Here, we split the pressure p 6 (see Fig. 2(b) ) into 



p £ =p £ +pI, 



(3.10) 
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with 



r i 



pKp) 



el{p(p*-S)+p'( P *-S)(p-p* + 5) 



if P < P* - 8, 



if p > p* — 5, 



(3.11) 



and 5 = . We treat the p £ component of the pressure force explicitly and the p\ 
component implicitly. The implicit treatment of the pressure force in (3.9) as well as that 
of the mass flux in (3.8) provides the AP property in the small Mach-number limit [28] . 
However, if the full pressure force p e is treated implicitly, some instabilities may appear 
|2"8] . Keeping the p £ part explicit removes these instabilities [22] • The explicit part p e 
is chosen to be p £ /2 except and in a neighborhood [p* — 5, p*} of p* with 5 adequately 
chosen. In this interval p e Q is a quadratic function. Matching of p E Q and of its derivatives 
up to second order is ensured at p* — 5. The choice 5 = e^- guarantees that Pq' (p* — 5) 
remains finite and consequently that p e itself remains finite over the whole range [0, p*} 
uniformly as e — > 0. 



Remark 3.2 When a background pressure Pb(p) prevails (see Remark 2.1), the splitting 
(3.10) is replaced by (2.21) i.e. we let p e (p) = Pb(p) and p\(p) =ep(p). 

Now, following [221 12H], we insert (3.9) into (3.8) to get an elliptic equation for p n+l l 2 . 

.,71+1/2 _ n n 

; " At\A x (pl(p n+1/2 )) 

(3.12) 



P" 



At 

V x ■ q n + AtcV 



+ At\A x (p £ (p n )). 



Instead of solving (3.12) for p 



n+l/2 



we solve it for p™ +1 ^ 2 and then, deduce p 



n+l/2 



by 



inverting the relation pi = p\(p), which is possible thanks to the monotonicity of p\ . This 
procedure is more stable and satisfies the congestion cons traint p n+l < p* automatically. 
Once p n+l / 2 is known, q n+l / 2 is computed by solving (3.9). 

We now consider the second (relaxation) step (3.6), (3.7). We first note from (3.6) 
that the density is unchanged: p n+1 = p n+l / 2 . Then, |f2(t)| 2 = \q(t)\ 2 /\p(t)\ 2 is a solution 
of: 



:dt\U\ 



'i - \n\ 2 )\fi\ 



which can be solved analytically. We easily deduce that: 



n+l 



n n+l/2 



Q n + l/2 



e p 



At 



fr+ 1/2 . 



(3.13) 
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When /5 — y 0, (3.13) leads to a simple renormalization of the velocity field: 

Q n+l/2 



n+1 



n n+l/2 



3.3 Fully discrete scheme 

The full space-time discretization is presented in the ID case. For the sake of completeness, 
the 2D discretization is given in Appendix 1. Since the second step does not involve any 
spatial derivative, we focus on the first compressible step. For simplicity, we consider the 
domain [0, 1] and a uniform spatial mesh of step Ax = jj, where M is a positive integer. 
We denote by Uj = (p™, g™) T the approximations of U = (p, q) T at time t n = nAt and 
positions Xj = jAx, for j — 0, 1, .. . , M. We discretize (3.8), (3.9) with a semi-implicit 



version of the local Lax-Friedrichs (or Rusanov) method [10] as follows: 



pT 12 



At 



1 

Ax 



Qj+i/2(Uj 
- Qj , 2 (E/ 



jjn Ti-n+1/2 TJ -n+l/2\ 



n+l/2 
1j 



At 



+ 



1 

Ax 
1 

2Ax~ 



[F j+1/2 (U?,U? +1 )- F^ 1/2 (Uy_ lt UV 



Aj>i(p 



n+l/2, 
J+l > 



0. 



where the fluxes are given by: 



Q 



i+1/2 



j-in 

^i+1/2 



1 

2 
1 
2 



n+l/2 



n+l/2 



2^+1/2^+1 ~~ Pj)> 



' Pi 



+ c 



(g?+l) : 
p?+l 



+ A^(p" +1 ) + Apg(p?) 



° j+i/2 Wi+i 9* 



(3.14) 



(3.15) 



(3.16) 
(3.17) 



The numerical fluxes consist of the sum of centered fluxes and decentering terms introduc- 
ing diffusion. The latter are evaluated explicitly while some of the former are evaluated 



implicitly. More precisely, the central part of the mass flux in (3.16) and of the p\ pres- 



sure contribution to the momentum flux in (3.15) are solved implicitly. By contrast, the 



central part of the velocity transport flux and of the Pq pressure flux in the momentum 



flux (3.17) are kept explicit. This level of implicitness is sufficient to guarantee the AP 
character of the scheme [2H] in the small Mach-number limit. The quantity CJ +1 ^ 2 is the 
local diffusion coefficient and is given by: 



C 



J+l/2 



max{q\qvi}, q 



p" 



+ W(c 2 -c) 



+ A(pg)'0tf?). (3.1* 
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It is defined as the local maximal characteristic speed related to the explicit pressure p £ . 
Therefore, it remains bounded as e goes to zero, thanks to the fact that p e Q itself remains 
uniformly bounded when e — > 0. By contrast, if a fully explicit scheme was used, C" in 
(3.18) would involve Po+p\ and would not remain bounded as e — > 0. The quantity C™ +1 / 2 



provides the numerical viscosity which is needed to ensure the stability of the scheme. 
Based on this discretization, we can apply the same strategy as described in section 



3.2 to get a discrete elliptic equation for p. We substitute (3.15) into (3.14) and obtain: 



n+l/2 

Pj 



At 



2Ax 



At 
AAx 2 



2Ax 
At 

2Ax 2 



[C i+ i/ 2 (p" +1 - p]) - Cj-y 2 (j>j - Pj-i)] 



j+3/2 



J+l/2 



1^-1/2 + F "-3/2] - 0. 



This equation is consistent with the time semi-discrete case (3.12). Then we get a non- 
linear equation for p\. 



p((Pi)j 



n+l/2^ 



At 2 



A 



' 8Ax 2 
P n ~ —(q n 



n+l/2 
i+2 



At 2 



Qi-i) + 



m 

At 



>r 1/2 + b e i)"-2 i/2 



4Ax 



[c i+V2 (^ +1 - # ) - c^/aC^ - 



+ 4 Ax 2 L i+3/2 J 



J+l/2 



1^-1/2 + ^3/2] 



We use Newton iterations to solve this nonlinear equation and get p™ +1//2 . The density 
pn+1/2 j g thgn obtained by inverting the nonlinear function p\ = p\ (p) with anoth er series 
of Newton iterations. Once p n+l l 2 is solved, q n+1 / 2 can be obtained by (3.15), which 
yields: 



q 



n+l/2 



with 



At 



A 



4Ax 
At 



n+l/2s 



fel/ 



J- 



2Ax L ^' +1/2 J - 1/2J 

Remark 3.3 (i) in the case where a background pressure prevails, we replace the splitting 



(3.10) by Pq(p) = Pb(p) and p\(p) = ep(p) (see Remark 3.2). 

(ii) In the case where c < 1, the previous scheme can be used as well. The only place 



where hyperbolicity is used is in the definition of the diffusion coefficient (3.18). In the 
case c < 1 and when the quantity inside one of the square roots is negative, we replace it 



by its absolute value, i.e. the second formula of (3.18) is replaced by 

1/2 

- U'-r) ~ + A(pg)'(p?) 



q 



p 



p 
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4 Numerical results 



4.1 Riemann problem test 



In the one-dimensional case, p and f2 only depend on a one-dimensional variable x and on 
time t. The vector field f2 can be described by its angle 9 with respect to the x-direction, 
i.e. f2 = (cos 9, sin 9). Then, the SOH system can be written as follows: 



p t + cos 9 p x — p sin 6* 9 X = 0, 
p £ '(p) 

9 t — A sin 9 p x + c cos 9 9 X 

P 



0. 



(4.1) 
(4.2) 



This system is hyperbolic with the characteristics speeds given by (2.15). However, this 



system is non-conservative, and shock waves are not uniquely defined [39]. A simple 
conservative form can be found 1271 and is written as follows: 



p t + (p cos 9 
d t fi(9)-d x 



x = 0, 
(of 2(9) 



Hp)) = o, 



(4.3) 
(4.4) 



where 



/i(0) = ln 



tan 



/ 2 (0) = ln|sin0|, g'(p) 



P £ '(p) 
P 



From this conservative formulation, the solutions of the Riemann problem can be analyti- 
cally derived. Indeed, by eliminating the shock speed in the Rankine-Hugoniot conditions 



associated to the conservative form (4.3), (4.4), the equation of the shock curves can be 
determined [ID]. These equations relate the right state (p r ,9 r ) to the left state (pi,9 e ) of 
the shocks. It is easy to see that for a given left state, there are two such curves and that 
they are given by the following equation: 



(p r - pe)(cf 2 (9 r ) - cf 2 {9i) - Xg(pr) + ^g(pe)) = (pcos(6> r ) - pecos(9 £ ))(fi(9 r ) 



This equation can be solved locally. 

However, conservative forms of the SOH model (4.1) 



A(9t)). 
(4.5) 



(4.2) are non unique. Other 
conservative forms than (4.3), (4.4) can possibly be found. Therefore, the analytical 



solutions of the Riemann problem derived from (4.5) may not be correct. 



m 



The SOH model provides an approximation of the Vicsek particle model, as shown 
[25] . Therefore, physically valid solutions of (4.1), (4.2) are those which are close 
to the underlying Vicsek particle model. In [43], it is shown that the relaxation model 
provides such a consistant approximation. By contrast, some other numerical methods 
(such as standard shock-capturing methods) select physically 'wrong' solutions, which 
are different from those of the Vicsek particle model. Therefore, the numerical method 
developed in the present paper, which relies on the relaxation model of [13], provides a 
consistant approximation of the physically correct solution. In order to determine whether 
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the analytic solutions of the Riemann problem based on (4.5) are the physically correct 



ones, we can compare them to a reference solution obtained by the present numerical 
method. 

We test a shock of the first family (i.e. associated to the smallest characteristic speed) 
under the condition c = A = 1. Since the conservative form requires that sin# ^ 0, the 
test example can be defined as follows: 

(p e , 9 e ) = (0.8, 0.14), (p r , r ) = (0.9969, 1.4502), 

with the shock speed deduced from the Rankine-Hugoniot relation and given by 



a 



-3.4136. 



(4.6) 



The numerical simulation also yields a shock wave, as can be seen on Fig. [4] obtained 
with parameters e = 10~ 4 , /3 = 10~ 7 , mesh sizes Ax = 0.005, At = 0.0005 at time 
t = 0.05. In these simulations, one-dimensional simulations have been performed using 
a two-dimensional code by imposing periodic boundary conditions along the top and 
bottom horizontal boundaries. The initial condition has been set independent of the y- 
coordinate. Fig. |4(a)| shows the density map in the two-dimensional domain, with a color 
coding (brown for p = 0.8 and black for p = 0.9969). Fig. 4(b) represents the vector 
field f2 with blue arrows. Arrows are almost horizontal for 9 = 0.14 and almost vertical 
for 9 = 1.4502. Fig. [5] displays cross-sections of the various components of the solution 
along the axis y = 0.5 as functions of x at times t — (left figures) and t = 0.5 (right 
figures) 



Figs. 5 (a 



and 5(b)| show the density p and figs. 5(c) and 5(d) show the x and 
y-components of the momentum q (top and bottom pictures respectively). By recording 
the position of the shock at time t = 0.5, it is possible to determine the physical shock 
speed. By inspection of Fig. |5j it is found to be a ~ —3.5714. This is different from the 
prediction (4.6) and the difference is significant, given the fine mesh size. Indeed, if the 
shock speed was given by (4.6), the shock would reach the boundary at time t = 0.1465. 



This is different from what we observe in Fig. [5] by more than 10 time steps. 

These numerical evidences show that the conservative form (4.3)-(4.4) does not provide 
the true physical solution, while the relaxation method, as shown in [43], does. 



4.2 2D simulations of collisions of herds 
4.2.1 Initial conditions 

In this section, we show simulation results in a truly 2D setting. Our main focus is the 
dynamics of the interaction of two clustered regions, which cannot be studied theoretically 
with the methods of [27], as explained at the end of section [2] The initial datum describes 
two clusters with densities close to the congestion density moving towards each other, 
inside a lower density fluid. Since the speed of the flow must be of norm one everywhere, 
the lower density fluid (background flow) is supposed to adopt a swirling motion around 
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t =0.05000 



t =0.05000 





(a) density p 



(b) velocity fl 



Figure 4: Shock-wave test problem. Density p (Fig. 4(a) ) and vector field f2 (Fig. 4(b) ) 
as functions of x and y at time t = 0.05. The density is color-coded, with color map 



indicated to the right of Fig. 4(a) The field f2 is indicated by blue arrows. 



0.2 0.4 0.6 0.8 



0.2 0.4 0.6 0.8 



(a) density p at t = 



(b) density p at t = 0.14 



1.4 0.6 0.8 

t =0.00000 



0.2 0.4 0.6 0.8 



(c) momentum q at t = 



0.2 0.4 0.6 0.8 1 

t =0.14000 



0.2 0.4 0.6 0.8 1 



(d) momentum q at t = 0.14 



Figure 5: Shock-wave test problem. Cross-section along the axis y = 0.5 of the density 
p (figs |5 (a)"] and [5(b) ) , s-component of the momentum q (figs 5(c) and 5(d), top) and y- 
component of momentum q (figs 5(c) and 5(d) , bottom) at time t = (left) and t — 0.14 
(right) as functions of x. 
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the center of the domain D = [0, 1] x [0, 1]. This leads to the following initial condition: 



n(aj,0) 

A -- 



p(x,0) = 0.8 x l AuB + 0.7 x 1 

D\(AUB), 

1\ . {-I s 



1 1 

6' 2 



1a 



x 



0/ 1 1b \ 



-D\( J 4uB) 



V(X 



1 2 

3' 3 



B 



1 10 
2' 12 



1 2 

3'3 



'-(</- o.sy 

v/(x - 0.5) 2 + {y - 0.5) 2 V ^-O- 5 '' 



where 1a denotes the indicator function of the set A. The definition of v corresponds 
to a background flow rotating around the center of the domain P = (0.5, 0.5)* counter- 
clockwise. The initial data can be seen in Figure [6] Fig 6(a) shows the density p in 



the two-dimensional domain, in color coding (larger densities are shown in darker brown 
color, and lower densities, in lighter brown color). The dark inner rectangle shows the 
large density clusters, surrounded by lower density regions in lighter brown. Fig. 6(b) 



provides a view of the vector field Q represented by blue arrows, in the two-dimensional 
domain. Inside the large density rectangle, the arrows are horizontal, and point towards 
the right in the left-hand half of the rectangle and to the left in the right-hand half of 
the rectangle. This shows that the inner rectangle is formed by two high density clusters 
moving one towards each other. The arrows outside the inner rectangle indicate that the 
flow is circulating counterclockwise around the inner rectangle. 
The numerical parameters are chosen as follows: 



10' 



P = 10' 



A = 1, Ax = 0.005, At = 0.0005. 



(4.7) 



The presence of a relatively large background density is necessary otherwise instabilities 
develop and generate negative pressures [22]. Developments are currently undertaken to 
prevent these instabilities and ensure the positivity of the pressure in all cases. We first 
investigate the case where c = 1, without background pressure (in this case, we use the 
decomposition (3.10) of the pressure). Then, we add a background pressure, and use 
the decomposition (2.21) instead. Then, we keep a background pressure, and investigate 
successively the cases c = 2 and c = 0.5. In the latter case 
not hyperbolic, and we use the method described in Remark 3.3 
appearance of complex valued characteristic speeds. 



the relaxation system is 
) to overcome the 



4.2.2 Case c = 1, without background pressure 

Figure [7] displays the density p and the vector field f2 at time t = 0.05, using the same 



representation as in figure As we can see on Fig 7(a) , the two clusters initially located 



in the inner rectangle are moving one towards each other and generate a congestion region 
(displayed in black color) between them. The congestion constraint generates a deflected 
wave which propagates in the vertical direction and compress the background flow, also 
bringing it close to the congestion density. This compression occurs in the central parts 
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of the top and bottom horizontal boundaries of the high density rectangle, as indicated 
by the black and dark brown colors in this areas. Finally, due to its rotation about 
the high density inner rectangle, the background flow impinges against the high density 
clusters along the top left and bottom right horizontal boundaries of the high density 
inner rectangle. This also produces some congestion there, as indicated by the black and 
dark brown colors in these areas. The two moving clusters leave a vacuum region behind 
them, which is shown by the light beige color along the outer vertical edges of the high 
density rectangle. Similarly, due to the rotation of the background flow, a vacuum region 
is created along the bottom left and top right horizontal boundaries of the high density 
inner rectangle shown in light beige color in these areas. 

Figure 7(b) highlights the geometrical features of the vector field Q. Except in the 
congested regions, the vector field is not affected. Therefore, the swirling pattern of the 
background flow away from the inner rectangle is clearly visible. Similarly, f2 is not 
changed in the lower left part of the left cluster and upper right part of the right cluster: 
it is directed in the horizontal direction and points to the right in the left cluster and 
to the left in the right cluster. However, f2 is strongly affected in the congested region. 
The flow abruptly changes direction towards a globally vertically directed flow in the 
central part of the clustered region. This flow is directed downwards in the left cluster 
and upwards in the right cluster. This is consistent with the direction of the swirling 
flow, which seems to force the circulation inside the congested region to adopt a similar 
counterclockwise rotation as for the background flow. As a result, an S-shaped shock line 
appears inside each cluster, separating the regions of unperturbed horizontal motion to 
those of vertical motion. The upwards and downwards vertical motions are separated by 
a slip line located in the middle between the two clusters. This slip line is almost vertical 
(slightly tilted in the south-east to north-west direction). This slip line ends up with 
two counterclockwise rotating vortices. These vortices are generated by the vertical flows 
inside the clusters that short-circuit the swirling motion. They are located where the 
vertical slip line rejoins the swirling motion at the upper left and lower right horizontal 
boundaries of the high density clusters. The vertical circulation which takes place in the 
inner part of the clusters penetrate the background flow on the sides of the vortices and 
generates a shock line with a hemi-circular shape. This hemi-circular shape reproduces 
the shape of the transition between congested and uncongested densities that can be 



oberved in Fig. 7(a) , near the central parts of the top and bottom horizontal boundaries 
of the high density inner rectangle. 

From these observations, we deduce that f2 is not smooth. The constraint = 1 
contributes to generating complex patterns. Their geometric features are outlined in 



section 2 and in Fig. 3(a) Indeed, as a consequence of the fact that Q is constant 



along straight lines normal to itself, its integral lines are parallel curves to each other. 



This feature is clearly seen on Fig. 7(b) However, at singularities, the parallel curves 



are interrupted (see Fig. 3(b)). Here, due to the collision of the two clusters, vortex 



singularities, discontinuities and slip lines form. The resulting topology takes a form 
that is reminicent of domain walls in micromagnetism. Due to the time dynamics, these 
domain walls evolve in the course of time and can change topology. 
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As time goes on, the two herds interact and form a larger congestion region that 
spreads towards the low density zones. This can be seen in Figure [8j which shows the 
solution at time t = 0.1 (see density p in Fig. 8(a) and vector field f2 in Fig. 8(b), with 
the same representation as in the previous figures). In Fig. 8(a)| we note that, due to 
wave dispersion, the densities inside the congested regions of the background flow and 



of the cluster tails have decreased below the congestion density (the color in Fig 8(a) 
is brown and not black in these regions). Congestion becomes strictly restricted to the 
compression region between the two clusters. The vacuum region following the tail of 
the clusters has also widened, because the clusters have moved further, leaving a region 



devoid of particles behind them. The general flow structure discussed above for Fig. 7(b) 
still remains on Fig 



Kb) 



However, the slip line at the interface between the two clusters 
in the middle of the high density region has been transformed into a saddle singularity. 
Simultaneously, the horizontal flow inside the clusters seem to be more strongly affected. A 
splay field structure appears near the center of the clusters close to the saddle singularity. 
There, the horizontal flow splits into upwards and downwards components. The roughly 
straight top and bottom horizontal boundaries of the congestion region, which can be 



seen in Fig. 8(a) as an abrupt transition from black to brown corresponds to a contact 
discontinuity. Indeed, in Fig. 



Kb) 



along the same line, there is no discontinuity of the 
field f2. Therefore, across this line, the density suffers an abrupt discontinuity while the 
velocity is continuous. Finally, the entire structure of the high density inner rectangle 
has rotated by an angle of about 10 degrees in the counterclockwise direction, i.e. in the 
same rotation direction as the background swirling flow. 



4.2.3 Case c = 1, with background pressure 

We now investigate how the presence of a background pressure modifies the density pat- 



terns and the flow structure. We change the pressure relation to (2.21) and we take pb 
as (2.22) with k = 1 (see also remark 2.1). The other parameters are given by (4.7). 



The simulation is more robust since the background pressure prevents the appearance of 
negative densities. 

The results at time t — 0.1 are shown in Fig. [9] (see density p in Fig. 9(a) and 
vector field f2 in Fig. 



9(b) 



with the same representation as in the previous figures). 
This figure must be compared with Fig. |8j where the results without background density 
are displayed, the other parameters being the same. We can see that the effect of the 
background pressure is to smear out the vacuum regions and to propagate the compression 
wave faster. Indeed, on Fig. 



9( a ; 



we remark that the vacuum regions have partially filled 
up due to the development of rarefaction waves at the outer vertical boundaries of the high 
density rectangle. These rarefaction waves bring the fluid from the backgound low density 
regions into the vacuum regions, as can be seen on Fig. 9(b)| Indeed, in these vacuum 
regions, the direction of the vector field f2 is oblique and points from the background 
flow towards the vacuum region. By contrast, the direction of f2 for the corresponding 
pressureless case in Fig. 8(b) is vertical. Like in the pressureless case (Fig. 8(a)), a 



vertical flow is created by the collision of the clusters. It arises from the deflection of the 
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cluster velocities into the vertical direction by the congestion constraint. This vertical flow 



impinges on the background flow and creates a hemi-circular wave- front (see Fig. 9(a) 



middle part of the top and bottom horizontal boundaries of the high-density rectangle). 
However, by contrast to the pressureless case, the wave front propagates further and has 
a more regular circular shape. This is obviously due to acoustic wave propagation in 
the low-density background flow, which conveys information about velocity and modifies 



the flow direction downstream. Indeed, in Fig. 9(b) , it can be seen that the direction 
of f2 is abruptly changed along the shock curve that bounds the region of compression. 
The acoustic wave conveys information about the vertical velocity of the cluster into the 
low density background. These changes of flow direction contribute to pre-compress the 
gas, far downstream the congestion region. Similar to the pressureless case, we notice 
that the nearly straight top and bottom horizontal boundaries of the congestion region 



(see transition from black to brown on Fig. 9(a)) are contact discontinuities. Indeed 



across this line, there appears no discontinuity of S7 on Fig. 9(b) The density in the 
cluster tails (Fig. |9(a)[ ) is smaller than in the pressureless case (Fig. 8(a)). This is 
also due to the development of rarefaction waves that contribute to fill in what was the 
vacuum region in the pressureless case. We notice that the two counterclockwise rotating 
vortices and the saddle singularities in the middle of the cluster interaction region remain 
(compare Figs. |9(b) and 8(b)). Therefore, the general structures of the flow are similar 
in the background or no-background pressure cases, but the details of these structures 
are quantitatively affected in a perceivable way by the presence or not of the background 
pressure. 



4.2.4 Case c ^ 1, with background pressure 

Influence of c: theoretical analysis. Before investigating the influence of the parame- 
ter c on the simulation results, we perform some analysis of the effect of this parameter on 
the propagation of waves. We first digress on the physical significance of the characteristic 



speeds (2.15). 



The characteristic speeds (2.15) are the propagation speeds of the waves generated 



by a small perturbation (5p,5Q) of a uniform state (p, SI), in the linear approximation. 
Since these waves are independent of the spatial coordinate normal to the propagation 



direction, it is legitimate to consider the one-dimensional system (4.1), (4.2). Introducing 
the variable u = cos# G [—1,1], this system can be written in terms of the unknowns 
(p, u) in matrix form as follows: 



A{p,u) 



with 



A{p,u) 



u 



-A(l 



u 



0. 



P 

CM 



(4.9) 



where x is the coordinate along the propagation direction and A = X p — . The quantity 
u is the projection of the vector f2 along the propagation direction. The characteristic 
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speeds £± are the two eigenvalues of A(p,u), i.e. the roots of the polynomial 



det(A(p,u) — £ Id). Their expression, deduced from (2.15), is recalled here for the sake of 
convenience: 

e± = ^ ((l + c)u±VA) , A = (1-c)V + 4A(1-m 2 ). (4.10) 

Now, the perturbation waves 5m) are not arbitrary, but are colinear to the asso- 
ciated eigenvectors. In other words, a simple wave propagating with velocity £± is such 
that 

5u 

where H± is the eigenvector of A associated to £±. This condition leads to the following 
relations between Sp and Su: 

-((1-c)utVA)<*P + p<5w = 0. (4.11) 



In (4.11), the minus sign is associated to £ + and vice-versa. 

The special case u = 1 is particularly simple. It corresponds to an unperturbed state 
such that the vector f2 is parallel to the propagation direction. In this case, we have 
£± = |((1 + c)u ± |c— To simplify further, we must discuss the position of c with 

respect to 1. 

- If c > 1, we have: 

6- = w, £+ = cu, 
and the associated waves are such that 

6u- = 0, p5u + = (c - 1) 5p+, (4.12) 

where the indices '+' and '— ' refer to the wave associated to £ + and 

- If c < 1, we have: 

f- = ctt, ^+ = u, 
and the associated waves are such that 

p 5 M _ = -(l- c )5p_, 5m+ = 0. (4.13) 

We notice that, in both cases, one characteristic speed coincides with the velocity u 
and that the associated wave is a contact discontinuity: the velocity is continuous on both 
sides of the waves and the density varies arbitrarily. 

The other characteristic speed corresponds to the velocity cu. Therefore, it is the 
velocity associated to the transport of f2. The corresponding eigenvector shows that 
velocity perturbations associated to this wave affect the density. Suppose that initially, 
5u is decreasing with respect to x. This corresponds to a situation where particles in the 
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rear are faster than particles in the front. In usual fluids, this contributes to stiffening the 
velocity front, ultimately leading to a shock, like in the Burgers equation. Simultaneously, 
a density bump forms in the region where the velocity fronts stiffens. Let us examine how 
density perturbations behave in the presence of a decreasing velocity perturbation profile 
in the situation where c ^ 1. In the case c > 1, the second relation (4.12) shows that 
5p + has the same sign as 5u + and is therefore, decreasing. Therefore, the influence of 
c > 1 is to lower the density increase due to the stiffening of the shock. In the case c > 1, 
we thus expect a lowering of the amount of congestion generated by the dynamics. By 
contrast, in the case c < 1, the first relation (4.13) shows that 5p- has the opposite sign 
to that of 5u- and is therefore, increasing. Therefore, the influence of c < 1 is to raise 
the density increase due to the stiffening of the shock. In the case c < 1, we thus expect 
an amplification of the amount of congestion generated by the dynamics. 

Case c > 1. We now illustrate this analysis by numerical simulations. In Fig. [TUJ 
we show simulation results obtained with c = 2, the other parameters being given by 
Again, we consider a background pressure given by (2.22) with k — 1. The 

with the same 



(4.7) 



density p is shown in Fig. 10 (a 



10(b) 



and the vector field f2, in Fig. 
representation as in the previous figures. This figure must be compared with Fig. [9j 
where the case c = 1 with the same background pressure is displayed. As expected from 
the analysis of the role played by the constant c above, the size of the congested area is 
smaller than in the case c = 1 (compare the surface occupied by the black regions on Fig. 



10(a) and Fig. 9(a)). The vertical flow created by the deflection of the clusters due to 
the congestion constraint impacts the background flow more deeply. The hemi-circular 
structure generated by this impact extends almost up to the top and bottom horizontal 
boundaries of the domain. But simultaneously, the area where the congestion density is 
reached is smaller. Therefore, most of the compression is due to the fast propagation of 
the velocity information, due to the large value c = 2. This propagation pre-compress the 
fluid downstream the congestion and contributes to reducing the amount of congestion. 
Again, we observe that the horizontal boundary of the congested region is a contact 
discontinuity as no associated discontinuity can be observed on the chart of f2, in Fig. 
10(b) The two counterclockwise vortices have transformed into unstable nodes as seen 
from the flow patterns in Fig. 10(b) , and simultaneously, the area near these unstable 



nodes has emptied from particles. This area appears in light beige color in Fig. 10(a)| 
indicating that the density is very small there. The central saddle point has remained 
and the congested area adopts a remarkable X-pattern around this point. The vacuum 
regions along the vertical boundaries and bottom left and top right horizontal boundaries 
of the initially high density inner rectangle have been almost completely filled and are 
hardly noticed on Fig. |10(a) 

Case c < 1. We now turn to the opposite case, namely the case c = 0.5. The case 
c < 1 is of broader significance than the case c > 1. The case c < 1 covers traffic systems 
and a larger class of biological swarming systems such as fish schools, mammal herds, 
. . . [251 IS]- By contrast, the case c > 1, which occurs with backwards vision (31], would 
apply e.g. to locust swarms. As discussed in Sec. |3.1[ the relaxation approximation is 
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11 



To bypass this problem, we adopt the 
we show the simulation results with 



not hyperbolic in general, although the SOH is 
strategy described in Remark 3.3 (ii). In Fig. 
c = 0.5 in the presence of a background pressure, the other parameters being given by 
(4.7). The density p is shown in Fig. 11(a) and the vector field fl, in Fig. [11(b) with the 
same representation as in the previous figures. This figure must be compared with Figs. [9] 
and 10 where the cases c = 1 and c = 2 with the same background pressure are displayed. 
The qualitative features of the solution in Figs. [9 10 and 11 are fairly similar. This is a 
strong indication that the solution obtained with c = 0.5 is correct and this validates the 
strategy consisting in using the relaxation system, in spite of its non-hyperbolic character. 
Unfortunately, there are no analytical solutions available for this model. So, we cannot 
provide a better validation than through this analogy for the time being. 

Again, the analysis of the role played by the constant c above is confirmed by the 

Indeed, the size of the congested areas is larger than in the 
we notice that there remains a congested area 



11 a 



11 a 



inspection of Fig. 
case c = 1. For instance, in Fig. 
(black color) in the part of the cluster which has not yet been deflected (region just below 
(resp. above) the upper (resp. lower) vortex), while the same region in Fig. 



9(a) 



is 

no more congested. The vertical flow created by the deflection of the clusters due to the 
congestion constraint penetrates the background flow less deeply. But, on the other hand, 
all the area included in the hemi-circular structure generated by this impact is very close 
to congestion (dark brown color on Fig. 11(a) ), while only one side of this region is close 
to congestion in Fig. 



9 a 



Propagation of velocity information is slow, due to the small 
value c = 0.5. Therefore, the pre-compression of the fluid downstream the congestion is 
less efficient. The uncompressed downstream fluid does not oppose to the propagation 
of the congestion and therefore, the amount of congestion is larger. Again, the feature 
that the horizontal boundary of the congested region is a contact discontinuity can still 
be observed by comparing the density chart (Fig. Fig. 11(a)) and the chart of f2 (Fig. 
11(b)). The two counterclockwise vortices are there (Fig. 10(b)), but the vacuum area 



close to these vortices is almost non-existing (see Fig. 11(a)). The central saddle point 
is also there. The vacuum regions along the vertical boundaries and bottom left and top 
right horizontal boundaries of the initially high density inner rectangle are fairly visible 
on Fig. 11(a) 



4.2.5 Conclusion of the two-dimensional tests 

From these tests, we first conclude on the validity of the numerical strategy developed 
here. This strategy, consisting of the combination of the relaxation model of |43j and 
of the AP-method of [22] proves its effectiveness, even in the cases where the relaxation 
model is not hyperbolic. The second conclusion of these tests is that the density and 
flow patterns generated by the SOH model can be very complex. The parameter c, which 
tunes the speed of propagation of velocity information with respect to the material flow 
propagation plays a key role in the type of structures that are generated. The lower value 
of c, the stiffer the dynamics is, with the generation of large congested areas and stiff 
density and velocity gradients. The constraint of the norm one velocity generates various 
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types of singularities such as vortices, nodes, slip lines, shock waves and saddle points. 
These structures are subject to topological changes as time evolves. It is also a remarkable 
feature that the method is able to reproduce these singularities and follow their evolution 
in time efficiently. 



t =0.00000 



t =0.00000 




0.2 0.4 0.6 




(a) density p 



(b) n 



Figure 6: Initial data for p (Fig. 6(a)) and f2 (Fig. 6(b)), as functions of x and y. The 



density is color-coded, with color map indicated to the right of Fig. 6(a) The field f2 is 
indicated by blue arrows. 



5 Application: a model of path formation in crowds 



In this section, we show how the relaxation model with congestion (3.1), (3.2) can be 
used to simulate path formation in crowds. To this aim, we extend it to two-fluid flows 
to model different groups of pedestrians heading towards opposite directions. We denote 
by (p+,q + ) (respectively (p_,qr_)) the density and momentum of pedestrians heading to 
the right (respectively to the left). The desired velocity of each group of pedestrians is 
denoted by w ± . The system satisfied by (p + , q + ,w + , p_, <7_,u>_) is written as follows: 

p +it + V x ■ q + = 0, (5.1) 

Q+,t + V, • + V,(P £ (P)) = ^(p+™+ - q+), (5.2) 

(p+w + ) t + V x ■ (w + ® q + ) = 0, (5.3) 
P-,t + V x -q_ = 0, (5.4) 

q , t + v, • ( q ~®_ q ~ ) + V,(p £ (p)) = - (5.5) 

(p_w_) t + V x - (w_®q_) =0, (5.6) 
P = P+ + P- (5-7) 
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Figure 7: Cluster collision. No background pressure and c = 1. Density p (Fig. 7(a) ) and 
velocity f2 (Fig. 7(b) ) at t = 0.05, as functions of x and y. The density is color-coded, 
with color map indicated to the right of Fig. 7(a) The field f2 is indicated by blue arrows. 





Figure 8: Cluster collision. No background pressure and c = 1. Density p (Fig. 8(a) ) and 
velocity f2 (Fig. 8(b) ) at t = 0.1, as functions of x and y. The density is color-coded, with 
color map indicated to the right of Fig. 8(a) The field f2 is indicated by blue arrows. 
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Figure 9: Cluster collision with background pressure and c = 1. Density p (Fig. 9(a) ) and 
velocity f2 (Fig. 9(b) ) at t = 0.1, as functions of x and y. The density is color-coded, with 
color map indicated to the right of Fig. 9(a) The field f2 is indicated by blue arrows. 




Figure 10: Cluster collision with background pressure and c = 2. Density p (Fig. 10(a) ) 
and velocity Q (Fig. 10(b)) at t = 0.1, as functions of x and y. The density is color- 
coded, with color map indicated to the right of Fig. 10(a) The field f2 is indicated by 
blue arrows. 
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t =0.10000 



t =0.10000 





Figure 11: Cluster collision with background pressure and c 
and velocity f2 (Fig 



0.5. Density p (Fig. |ll(a)D 
11(b)) at t = 0.1, as functions of x and y. The density is color- 

The field Q is indicated by 



coded, with color map indicated to the right of Fig. 
blue arrows. 



11(a) 



Eqs. (5.1 ), (5.4) are the continuity equations and eqs. (5.2), (5.5), the momentum balance 
eqs. for each species of pedestrians. The left-hand side of these equations form two systems 
similar to the compressible isentropic gas dynamics equations. However, they are coupled 



by a single pressure p £ (p) inside the momentum balance equation, where p given by (5.7) 
is the total density of pedestrians. The pressure has a stiff profile when the density reaches 
the congestion density p* given by (2.4), (2.5). The tendency for pedestrians to rejoin 



their target velocity is modeled by the relaxation terms at the right-hand sides of the 
momentum balance eqs. (5.2), (5.5), with relaxation rate /3 . The desired velocity is a 



can be written in non-conservative forms as 



0. 



u± 



quantity attached to the pedestrians, and is therefore passively transported by the flow. 
This is expressed by eqs. (5.3), (5.6). Indeed, using the continuity eqs. (5.1), (5.4), these 
eqs. 

s± 

which express a passive transport of w± by the velocity u±. This model bears analogies 
with the the Aw-Rascle model of car traffic [31 18] and its extension to pedestrian traffic 

m. 

The following test problem shows how two flows with opposite desired directions in- 
teract. The initial data correspond to a flow at rest with a random perturbation on the 
density. The steady flow at rest is defined as follows: 



Ps 



0.4, 



\x 



0. 



w. 



x 



[X 



0.4, q s _(x) = 0, w s _(x) 



where pedestrians heading to the right (resp. left) have desired velocity along the x axis 
pointing in the positive (resp. negative) direction. We perturb the density in the square 
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[1/3,2/3] x [1/3,2/3] by a random noise to model some inhomogeneities in the crowd. 
The perturbation p T+ is sampled out of a uniform distribution in the interval [—0.19, 0.19] 
and is constant on squares grouping 25 cells of the simulation mesh (i.e. squares made of 
5 spatial steps on each side). Finally, the initial conditions are set up such that 

p+(x,0) = Ps+{x) + p r +(x), p(x,0) = 0.8, p_(x, 0) = p(x, 0) - p+(x, 0), 
q±(x,Q) = q s± (x), w±(x, 0) = w a ±(x). 

The non-uniform initial densities allow for some non-uniform motion to develop. The 
low density regions give room to pedestrians to pass through while, due to the den- 
sity constraint, the high density regions act as bottlenecks preventing pedestrians to go 
through. To illustrate the process more clearly, we choose a relatively slow relaxation 
process /3 = 0.5. The mesh size is chosen to be such that Ax = Ay = 0.005 and the 
time step is At = 0.0005. Fig. 12 provides a representation of the initial densities of 
the right-going pedestrians p + (fig. 12(a)) and of the left-going ones p_ (fig. 12(b)) in 



the two-dimensional domain, using a color code. The perturbation of the uniform initial 
densities inside the inner rectangle takes the form of a checkerboard with cells of random 
color ranging from light beige color (small densities) to dark brown ones (large densities). 

To solve this system, the splitting method is used in the same way as described in 
section [3] During the conservative step, the equations for (p, q) and w are decoupled and 
are written as follows: 



p+,t + Va; • q + = 0, 

Q+ ® q. 



t+,t 



v.(p6(p)) + v.(p5(/0) = o, 



(p+w + ) t + V x ■ (w + <g> q + ) = 
P-,t + Vcc ■ q = 0, 
q ® q 



Q ,t + Va, • 



P- 



v,(p^(p)) + v 3 ,(p £ 1 (p)) = o, 



{p_w_) t + V x - (w^®q_) =0, 



(5.8) 

(5.9) 

(5.10) 
(5.11) 

(5.12) 

(5.13) 



where the pressure is decomposed into p £ (p) = p%(p) +pf(p) as in (3.10), (3.11). We can 



add (5.8) to (5.11) and (5.9) to (5.12) and get the following conservation eqs. for the total 
density p and total momentum q = q + + q as follows: 



Pt + • q = 0, 



, 2V a Mp))+2V a! (K(p))=0. 

P+ J V P- 

In a similar way as in the previous sections, we can implement the AP schemes for this 
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system. The semi-discretization can be written as follows: 



p n+l _ p n 

At 
q n+l _ q n 

At 



+ V x q 



n+l 



+ v 3 



(5.14) 

+ 2V x {pliP n )) 
+2V x (pl(p n+l ))=0. (5.15) 



The s trateg y for s olving this system is the same as before. The elimination of q n+l between 
(5.14) and (5.15) leads to an elliptic equation for the pressure p\' n+ ■ From Pi' n+ , the 
total density p n+1 can be obtained by inverting the relation p\(p n+l ) 



p{ . 



Then, the 



momenta q± can be updated thanks to the discrete, semi-implicit versions of (5.9) and 



(5.12). Finally, the densities p± can be updated thanks to the discrete, implicit version of 



(5.8) and (5.11). The equations for w± (5.10) and (5.13) are decoupled and can be solved 



easily by standard shock capturing methods. We leave the details to the reader. 

The second step of the splitting method is the relaxation step. It is written as follows: 

1 . 

p±,t = 0, w± >t = 0, q± tt = -n(p±w±- q±), 

and can be solved easily. 

The right-going pedestrian density p + is shown in Fig. 13 at two different times (fig. 
t = 0.025 ; fig. 13(b) t = 0.05) in the two-dimensional domain, using a color 
Due to the relaxation to the desired velocities, the flow is set into motion. The 

see Fig. |l2(aj| 



code. 

sharp edges between the checkerboard squares in the initial density p + 
progressively fade away in Figs. 13(a) and |13(b) as the two species of pedestrians mix 
together. 

In order to highlight the flow structure, we compute the following two quantities which 
measure which of the left or right going flows is dominant: 



Dp = p A 



P- 



where q± x is the first component of q±. The quantity Dp measures which pedestrian 
species is dominant (Dp > if there are more right-going than left-going pedestrians at 
this point, and vice- versa if Dp < 0). The quantity Dq\ measures which flow is dominant 
[Dqi > if the right-going flow is dominant and vice- versa if Dqi < 0). We expect 
that Dp and Dqi fluctuate about zero since initially p + ~ p_ up to random fluctuations 
and the initial desired velocities are opposite. But the morphology of the fluctuations 
results from the dynamics. In particular, we expect that alternating regions where right- 
or left-going pedestrians dominate will form. In these regions, the dominant flow will be 
right-going or left-going respectively. 

In Fig. 14 and 15, Dp and Dqi are plotted at different times as functions of x and 





y using a color code (Fig. 14 (a 
Fig. [15(b)] t = 0.075). 



t 



Fig. 14(b) t = 0.025 ; Fig. 15(a 



t = 0.05 ; 

In each of these figures, the top picture shows Dqi and the 
bottom picture, Dp. The color code is indicated by the color bar next to the figures 
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(the zero value is green, red colors show positive values, while blue colors, negative ones). 
This figures show the emergence of a phenomenon of path formation. Paths appear as 
elongated patches of the same color in Dqi (top pictures). Blue (resp. red) regions 
materialize paths of pedestrians moving from right to left (resp. from left to right). Path 
formation indeed results from the balanced influence of the relaxation towards the target 
velocity which forces pedestrians to move in one given direction and congestion avoidance 
which drives pedestrians towards the regions of lower density. Alternating regions of high 
densities of right-going and left-going pedestrians are materialized in the bottom figures 
by patches of red and blue colors. These patches reproduce the structure of the patches 
of Dqi but not exactly, meaning that some pedestrians are forced to adopt a motion in 
the transverse direction (y-direction) because their desired motion is impeded by large 
density clusters of opposite pedestrians. These simulations show the ability of the model 
to generate paths from initial random density fluctuations. 





(a) Density p + at t = 



(b) Density at t = 



Figure 12: Crowd model: Initial densities p + (Fig. 12(a)), p_ (Fig. 12(b)) as functions 



of x and y. The densities are color-coded, with color map indicated to the right of the 
figures. 



6 Conclusion 

In this paper, we have investigated the Self-Organized Hydrodynamic (SOH) model de- 
rived in [25J. Short-range repulsion is modeled by a singular pressure which becomes 
infinite at the congestion density. The singular limit of an infinite pressure stiffness leads 
to phase transitions from compressible to incompressible dynamics. We have proposed 
an Asymptotic-Preserving scheme which allows to treat the singular pressure with a very 
small but finite value of the stiffness parameter. We have performed numerical simula- 
tions which illustrate the efficiency of the scheme to treat the occurrence of congestions. 
A two-fluid variant has been proposed to model path formation in crowds. 

Future work will be devoted to the physical enrichment of the model. It will make it 
quantitatively more accurate for the modeling of herds of gregarious animals like ovine 



31 



t =0.02500 



t =0.05000 




0.475 




0.39 



(a) density p + at t = 0.025 



(b) density p+ at t = 0.05 



Figure 13: Crowd model: Density function of x and y at times t = 0.025 (Fig. 

l"3(a)T ) and t = 0.05 (Fig. [13(b)] ). The densities are color-coded, with color map indicated 
to the right of the figures. 
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(a) Dql and Dp at t = 



(b) Dql and Dp at i = 0.025 



Figure 14: Crowd model: (a)-(b) x-component of the momentum difference Dq\ (top) and 



density difference Dp (bottom) as functions of x and y at times t — (Fig. 14(a) ) and 



i = 0.025 (Fig. 14(b)). The numerical values are color-coded, with color map indicated 



to the right of the figures. 
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t =0.05000 
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(a) Dql and Dp at t = 0.05 



(b) Dql and Dp at t = 0.075 



Figure 15: Crowd model: (a)-(b) x-component of the momentum difference Dqi (top) and 



density difference Dp (bottom) as functions of x and y at times t = 0.05 (Fig. 15(a) ) and 



t = 0.075 (Fig. 15(b)). The numerical values are color-coded, with color map indicated 



to the right of the figures. 
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or bovine. It will also allow a more precise description of pedestrian flow and crowd 
behavior. To this aim, more detailed information from ethology and cognitive sciences 
will be included in the model. Another direction of improvement concerns the ability 
of the numerical method to treat vacuum regions and to provide sharp transitions from 
uncongested to congested areas free of any spurious oscillations. 



Appendix 1: The two dimensional full time and space 
discretization 

We consider the 2D case with domain Q = [0,1] x [0,1]. We denote by (xi,yj) = 
(iAxJAy), i = 0, • • • ,M X , j = 0, ••• ,M 2 , where M x = 1/Ax, M 2 = 1/Ay. Let 
U = (p,q) T , q = (qi,q2) T and Uij = U(xj,yj). To simplify the exposition, we define 



nu) = \ c * + ±° {p) ), g(u)-I 



Then the left-hand side of (3.4) and (3.5) can be written as 

dtp + d x qi + d y q 2 = 0, 

d t q + d x F(U) + d y G(U) + XV x (pl(p)) = 0. 

We denote by 



p 

e\n+l 



DIM{P)T +1 \DlM(p n+1 )) ' 



where Df^u, D\jU are the centered difference operators, which, for any scalar function u, 
are defined as follows: 

n x _ yn+hj ~ u i-l,j n y _ u i,j+l ~ u i,j-l 
l ' jU ~ 2Ax ' ^ 2Ay 

We also define the eigenvalues of the Jacobian matrix for two one-dimensional hyperbolic 
system as follows: 



2 / „2 



A« = c * ± J{# - c)% + A(p6)'(p), X {2) = c*,c* ± J(<* - c)% + A(pg)'(p). 

p p V p p p V p 

With these notations, the full discretization of the scheme takes the following form in 2D: 



A* Ai V i+ 2-J ^i-sd) Ay 



2 " 2 



A* Aa; V l+ 5 J % ~2^ J Ay\ 
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where the fluxes are 



jpn 



G n 



i,3 + i 



and 



^i+sJ = max {l A £' I' l\-+ijl> I' l\-+ijl}> 



( 2 ) I l\(2) 



C7 y+ i =max{|\W|J\W + J,|AjJ|,|AjJ +1 |}. 

Similarly to the one-dimensional case, by inserting the momentum equation into the 
density equation, we can get the following discrete elliptic equation: 



(2), , x (2) 



o n+l 



At 2 



1 



4 1 Ax 2 
1 



+^{^ [(^ + s/ 2 ,) (i) - (^ +1/2 .) (i) - {Fu„,r + {Fu^n 

1 



+ 



Ax Ay 



)(2) 



J-1/2J+1 



)(2) 



-l) (2) 



i-l/2j- 



-l) (2) ] 



+ 
+ 



+ ^ [(^ + 3/2) (2) - (^- +1/2 ) (2) - (G?„-_ 1/2 )« + (G?„-_ s/2 ) (2) ] 
,j(pi+l,j - p r i,j) - c i~l,j(pi,j - Pi-lj) 



2Ax 
At 
2Ay 



Here (-F?+i/2j)^ * s the fi rs ^ component of the vector F™ +1 i 2 ,j- Also, like in the one- 
dimensional case, we solve the above elliptic equation to get first then p n+l . Once 
pH+i j g known, we can g e t q n+l explicitly by solving 



q n+1 



At 



Ax\ J Ay 



At 



( G l3 + l- G l3-^-^*MY +1 - 
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